library(brms)
library(dplyr)
library(ggplot2)
library(ggthemes)
library(kfigr)
library(knitr)
library(patchwork)

The data set was read in as is, the columns subj, arm, cond, series and no were converted to factors, the column filename was deleted. A new column was added to the dataset, patient, based on the subj prefix (‘P’ or ‘S’), with a ‘0’ coding for healthy controls, and a ‘1’ for patients.

read.csv("features.csv") %>%
  mutate(filename = NULL,
         patient = factor(if_else(grepl("^P", 
                                        subj), 
                                  1L, 
                                  0L)),
         series = factor(series),
         subj = factor(subj),
         arm = factor(arm),
         cond = factor(cond),
         no = factor(no)) ->
  drum_beats

1 Exploration

1.1 Instruction Condition, and Arm Used

In a first attempt to get an overview, the descriptors’ density estimates were plotted, broken down by arm and cond.

dep_vars <- names(drum_beats)[which(!(names(drum_beats) %in% c("subj", "arm", "cond", "no", "series", "patient")))]
n_vars <- length(dep_vars)
plot_lst <- vector("list", length = n_vars)
plt_cnt <- 1
for (dv in dep_vars) {
  if (plt_cnt == n_vars) {
    # plt_pos <- c(2, 0.5)
    plt_pos <- "right"
  } else {
    plt_pos <- "none"
  }
  if (plt_cnt %% 4 == 1) {
    ylab <- "density"
  } else {
    ylab <- ""
  }
  plot_lst[[plt_cnt]] <- ggplot(drum_beats, 
                                aes(x = .data[[dv]],
                                color = .data[["cond"]])) +
    geom_density(na.rm = TRUE) + 
    scale_color_colorblind() +
    ylab(ylab) +
    geom_rug(alpha = 0.25) + 
    theme(legend.position = plt_pos) +
    facet_wrap(~ arm)
  plt_cnt <- plt_cnt + 1
}
print(wrap_plots(plot_lst, 
                 ncol = 4) +
        plot_annotation(
          tag_levels = "A"))

Fig. 1. Descriptor densities along with raw data points (rug ticks on x axes), broken down by arm (dominant [D] vs non-dominant [ND]) and instruction condition (controlled = C, normal = N).


Several peculiarities are immediately apparent from the plots in Fig. 1:

  • There are no substantial differences in the density distributions between normal (N) and controlled (C) strokes
  • attDur, LAT, and attFlat have very few unique data values (see rug ticks at bottom of plots)
  • Looking at the dominant arm alone (and excluding the vars attDur, LAT, and attFlat for the reason just mentioned), the densities have a trend toward a larger variance in the controlled condition (in other words: a less pointed distribution) with additional or larger humps in the tails, respectively

The lack of any clear-cut differences in the raw data between conditions suggests that it will be hard to find any differences by modeling.

The limited number of unique values in the attack-related measures (e.g. attDur only has 35 unique values in a total of 1102 observations across all subjects, conditions, and sides; that’s only 3 percent!) suggests that rounding errors propagated through the calculations. This is probably due to very similar, i.e. highly automated, and short attack times combined with the given sampling frequency, resulting in few data points, which in course of the calculations result in the observed phenomenon. But I’m just guessing here, as I do not know anything about these descriptors and how they are calculated or interpreted. Judging by the range of values of, i.e., attDur (2.59, 9.8ms, across both sides and conditions) and a sampling frequency of 48 kHz, this leaves us at 48 * (9.8 - 2.59) = 346 points to choose beginning and end of the attack. Given the supposed highly automatized motor program used to initiate the stroke, along with the laws of physics at play here (no pun intended), it is not very surprising to see very few unique values. My limited knowledge—or rather, my ignorance of anything sound-related—aside, from a statistical standpoint these measures do not seem suited to describe any differences between the experimental conditions investigated here.

Given the very few data points in the attack phase, any descriptor derived from such a short period of time (.e.g.. attSPL) cannot be judged as being stable in the sense of being reproducible. Hence I suggest to drop all attack-derived descriptors.

# drum_beats %>%
#   mutate(attDur = NULL,
#          LAT = NULL,
#          attSPL = NULL,
#          attSC = NULL,
#          attFlat = NULL,
#          attSpecFlat = NULL) ->
#   drum_beats

The increase in variance for the conditions N < C does not come as a surprise as I assume normal also means highly trained and thus automatized, whereas controlled involves less automatization and more ‘individualness’ both within and between subjects.

1.1.1 Conclusions

Looking at the above investigated descriptors and Danielsen et al. (2015), it seems reasonable to limit the modeling attempts to totDur, totSPL, totSC, and TC.

But Sofia wrote on 2020-07-09: “Francesco and I have discussed a bit related to descriptors and we want to concentrate on the “transient” period (although the name probably will change). Spectral Centroid (transSC) should be one, and I suggest transFlat for the other. Francesco, does that sound reasonable?”

On 2020-07-22 both Sofia and Francesco agreed upon transSC, transFlat, and transCrest as the probably most important response variables to look at.

1.1.2 Francesco’s comments

Update:

  • The lack of data points for attFlat is fixed by solving a bug in the feature extraction: now the downsampling ratio for the envelope extraction is reduced, and makes the MIRToolbox algorithm less sensitive to frames with an rms value of 0 (which returns an incorrect value of 0, since we have a geometric mean at the numerator). Now we have a larger spread which allows high values (1 = peaky envelope, 0 = smooth/flat envelope).
  • The crest factor is now calculated for the separate phases.

Regarding the attack phase: I agree with the remarks regarding duration. Given the fact that we are not comparing different instruments, I wouldn’t have expected a large variation. This is not surprising if we consider that our system is changing only slightly (same rototom, same drumstick, same action, a bit different tuning across subjects): in fact, the perceived timbral differences are so small that we are in trouble guessing on the descriptors.

The sampling frequency is even lower (\(f_s = 44100\) Hz). Even if we had a higher sample rate which could reveal some discrepancies in the attack durations, we would have to prove that they are perceptually relevant.

Therefore, I am happy to discard attDur and LAT (which is obviously a log-transformed duration, only there for the sake of consistency with the literature).

I am a bit more in doubt when it comes to discarding all the attack descriptors. Even if what Michael says is true from a statistical point of view, we should still be able to discriminate timbre on short time windows due to the high temporal resolution of our hearing. Attack phase descriptors (with the same definition of attack that we are using, which is most likely not coincident with perceptual attack) are employed in Câmara et al. (2020), and the Oslo group has a paper under construction which analyzes drum sounds in a similar manner (see OsloPlots.png in the OneDrive folder).

I am worried that merely taking the overall descriptors into account would introduce a lot of unnecessary and perceptually catching information — mostly the tonal part of the signal, i.e. the drum ringing in the last part of the decay. That’s why Sofia and I are suggesting to look at what we could call “main energy” or “early decay” phase (i.e. from max peak to temporal centroid).

Would it be feasible to set up 4 different models (i.e. one for each phase), at least in the univariate version?

As for the descriptors to include: although TC is employed in Danielsen et al. (2015), the PDFs are even more similar. I would go for Dur (except attack?), SPL, SC, and one between Flat, SpecFlat or Crest.

My (informal and biased) listening tests tell me that, at least for some subjects, I hear a pattern going towards a harsher (controlled) vs smoother (normal) timbre exactly at the hit point, plus slightly less (controlled) or more pitch/amplitude fluctuation. Hopefully this could be catched by spectral centroid, specrtral/temporal flatness, or crest factor. This should be independent of SPL unless the subject misinterpreted the instructions, therefore SPL acts as a sort of control variable in our model.

1.1.3 Model Considerations

There are several points that need to be considered before making a decision regarding the type of modeling to be done in this study, (1) the sample size places restrictions on the external validity; (2) data from small samples can be better modeled when regularization is in place to ‘tame’ the estimates; (3) the hierarchical structure of the data (subjects played several trials with either their dominant or non-dominant arm under two conditions) suggests a multilevel analysis of the data which would, in addition to the Bayesian regularization via priors, also results in shrinkage of the estimates; (4) given the small distribution differences between the two experimental conditions, along with the

  1. Given the extremely small sample size, any attempt to fit the data using traditional statistics, a.k.a. null-hypothesis significance testing, would make us vulnerable to all kinds of criticism. I therefore suggest to use Bayesian statistics instead, as it allows us to estimate probability distributions of parameters rather than confidence intervals around point estimates, and thus embraces uncertainty in estimates.
  2. Additionally, Bayesian regression uses prior probability distributions to arrive at sensible estimates. These priors regularize estimates, or draw them toward zero, a desirable effect which has long been recognized even in traditional statistics (e.g. ridge regression and lasso; Tibshirani (1996)).
  3. ANOVA (or regression with grouping variables) requires each subject’s trials to be averaged, say, within conditions or groups (or both), to be able to assess treatment effects; this results in the individual variation within a subject to get lost while it might have been very informative to include it in the analysis. Multilevel (or hierarchical) modeling allows to include the entire data structure (Fig. 2) so that no information gets lost through averaging, but all variation (both individual and treatment-driven) propagates through the analysis and ends up in the final results, allowing for more realistic credibility margins around estimates.
  4. Variance in traditional statistics is considered to be fixed (think homoscedasticity in ANOVA, or the \(\epsilon\) in regression formulas such as \(y_i = \beta_0 + \beta_1x_i + \epsilon\)), whereas it is an estimated quantity in Bayesian statistics and therefore can vary between conditions (or groups or subjects etc.) and thus capture varying variance in different conditions, or groups etc.

# chunk intentionally empty  
Group:                    Patient                       Healthy Control
                        /         \                     /               \
Instruction:    Normal            Controlled         Normal         Controlled
                | |    \            / |    \
Player:         1 2 ... n          1  2 ... n
               / \
Side    dominant non-dominant
        / | \       / | \
Trial  1 ..  p     1  .. p

Fig. 2. The multilevel structure of a data set should be reflected in the analysis.


1.1.4 Conclusions

Will start with a simple univariate model, add predictors and interactions, then a bivariate model, and finally a quadruple-variate model and see where this leads us.

1.2 Individual Responses

The humps in the density distributions in section Instruction Condition, and Arm Used made me curious where they might originate from. So in the following graph the density estimates of one of the descriptors, transSC, are plotted broken down by subj. and then also by cond, separately for patients and healthy subjects.

ggplot(drum_beats, aes(transSC, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() 

Fig. 3. transSC density distribution of participants, broken down by instruction condition.


It is obvious that some participants do not differ substantially between the normal and the controlled condition, whereas others do, and even markedly so. Additionally, there seems to be quite a spread of the centers of distributions across a wide range of values, suggesting very individual drum sounds.

The wide distribution of centers of mass between individuals made me want to further break down the plot.

ggplot(drum_beats, aes(transSC, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ patient)

Fig. 4. transSC density distribution of participants, broken down by instruction condition and group (patients and healthy subjects).


The above plot is interesting in that it seems to show that the four healthy subjects were more uniform than the patients in their transSC distributions and also, with the exception of S2, had very similar transSC distributions for the normal and the controlled conditions. In the patients, two had very similar distributions in both conditions (P1 and P5), whereas the two others showed differing results for the two conditions. So with regard to Sofia’s hypothesis (ch. Modeling) I’d argue, at least for transSC alone, having drummers play normal and controlled strokes would not allow subjects to tell the difference in a listening test. But maybe it would qualify as a screening test for movement disorders in drummers. Just a random thought.

While breaking it down it occurred to me that looking at individual variation (that is, between series) might also be enlightening.

ggplot(drum_beats, aes(transSC, 
                       color = cond, 
                       group = series)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ subj, nrow = 2)

Fig. 5. Density plots of individual series.


And it was! The plot in Fig. 5 made visible that P1 and P5 had consistently very low transSC values. P3 was consistent within the normal instruction condition, but had, on average, higher values, and with a lot more variation, in the controlled condition. P4 had more variation in both conditions, and higher transSC values in both; in other words, P4 was consistantly bad. (But then again: what do I know what bad is wrt transSC!)

The healthy subjects showed comparable variation and centers of mass within conditions, and all but one (S2) also across conditions.

Comparing the two rows of panels in Fig. 5 reveals that healthy subjects have more variation than the patients.

Although tempting, we probably shouldn’t get carried away and generalize to the populations of healthy drummers and ones with movement disorders, respectively.

For completeness’ sake the same exploratory analysis was carried out for transFlat and transCrest`:

ggplot(drum_beats, aes(transFlat, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() 

Fig. 6. transFlat density distribution of participants, broken down by instruction condition.


ggplot(drum_beats, aes(transFlat, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ patient)

Fig. 7. transFlat density distribution of participants, broken down by instruction condition and group (patients and healthy subjects).


ggplot(drum_beats, aes(transFlat, 
                       color = cond, 
                       group = series)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ subj, nrow = 2)

Fig. 8. transFlat density plots of individual series.


ggplot(drum_beats, aes(transCrest, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() 

Fig. 9. transCrest density distribution of participants, broken down by instruction condition.


ggplot(drum_beats, aes(transCrest, 
                       color = subj, 
                       linetype = cond)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ patient)

Fig. 10. transCrest density distribution of participants, broken down by instruction condition and group (patients and healthy subjects).


ggplot(drum_beats, aes(transCrest, 
                       color = cond, 
                       group = series)) + 
  geom_density(alpha = 0.1) +
  scale_color_colorblind() +
  facet_wrap(~ subj, nrow = 2)

Fig. 11. transCrest density plots of individual series.


1.3 Covariance Matrix

new_names <- drum_beats
names(new_names) <- gsub(pattern = "trans", 
                         replacement = "earlyDecay", 
                         x = names(new_names), 
                         fixed = TRUE)
cov_mat <- cov(new_names[-c(1:5, ncol(new_names))], 
               use = "pairwise.complete.obs")
kable(cov_mat, digits = 2)
totDur attDur decDur earlyDecayDur LAT totSPL attSPL decSPL earlyDecaySPL totSC attSC decSC earlyDecaySC TC totFlat attFlat decFlat earlyDecayFlat totSpecFlat attSpecFlat decSpecFlat earlyDecaySpecFlat totCrest attCrest decCrest earlyDecayCrest
totDur 8580.54 32.18 7641.34 907.01 5.93 -125.00 -205.08 -41.76 -90.45 -449.09 -7391.47 1448.04 -6574.23 999.35 -0.79 1.26 -2.41 2.57 1.05 -0.17 1.31 0.12 -90.78 -13.10 27.04 -39.39
attDur 32.18 1.15 30.22 0.81 0.23 0.26 -0.66 0.25 0.43 8.36 -171.34 47.75 -101.81 2.68 -0.02 0.02 -0.02 0.01 0.01 0.00 0.01 0.00 -1.26 -0.26 -0.03 -0.40
decDur 7641.34 30.22 6838.06 773.06 5.66 -108.21 -177.56 -36.46 -76.60 -325.99 -7202.60 1404.24 -5920.35 867.57 -0.87 1.19 -2.28 2.27 0.94 -0.17 1.16 0.11 -82.96 -12.27 25.25 -35.30
earlyDecayDur 907.01 0.81 773.06 133.13 0.04 -17.05 -26.86 -5.55 -14.28 -131.46 -17.53 -3.96 -552.06 129.09 0.09 0.05 -0.10 0.28 0.11 0.01 0.13 0.02 -6.57 -0.57 1.83 -3.68
LAT 5.93 0.23 5.66 0.04 0.05 0.10 -0.08 0.10 0.14 1.11 -37.72 8.53 -19.80 0.40 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 0.00 -0.26 -0.05 -0.01 -0.08
totSPL -125.00 0.26 -108.21 -17.05 0.10 35.41 35.05 34.35 35.05 -118.69 -128.44 -209.89 -29.55 -18.56 0.03 0.01 0.07 -0.02 -0.03 0.00 -0.04 -0.01 -0.25 0.01 -0.77 0.40
attSPL -205.08 -0.66 -177.56 -26.86 -0.08 35.05 37.49 33.06 34.36 -119.69 94.87 -274.74 122.55 -30.48 0.04 -0.03 0.10 -0.06 -0.05 0.00 -0.06 -0.01 2.08 0.38 -0.61 1.09
decSPL -41.76 0.25 -36.46 -5.55 0.10 34.35 33.06 34.70 34.18 -130.19 -122.66 -211.16 -48.42 -5.80 0.03 0.02 0.04 0.01 -0.02 0.00 -0.03 0.00 -1.13 -0.10 -0.53 0.06
earlyDecaySPL -90.45 0.43 -76.60 -14.28 0.14 35.05 34.36 34.18 34.88 -119.50 -177.76 -202.01 -60.17 -15.25 0.02 0.02 0.06 -0.01 -0.02 -0.01 -0.03 -0.01 -0.64 -0.06 -0.63 0.24
totSC -449.09 8.36 -325.99 -131.46 1.11 -118.69 -119.69 -130.19 -119.50 5907.01 1707.28 8098.86 1471.12 -94.91 -0.64 0.24 -0.75 -0.58 0.58 0.17 0.72 0.12 9.18 -3.27 1.29 5.57
attSC -7391.47 -171.34 -7202.60 -17.53 -37.72 -128.44 94.87 -122.66 -177.76 1707.28 62304.01 -5496.05 21851.87 -301.21 5.33 -5.54 6.07 -4.17 -1.74 1.92 -2.23 -0.10 314.17 49.56 -7.49 83.34
decSC 1448.04 47.75 1404.24 -3.96 8.53 -209.89 -274.74 -211.16 -202.01 8098.86 -5496.05 13980.47 -3473.23 115.29 -1.76 1.22 -2.51 0.47 1.42 0.13 1.80 0.24 -53.07 -16.77 -0.43 -13.80
earlyDecaySC -6574.23 -101.81 -5920.35 -552.06 -19.80 -29.55 122.55 -48.42 -60.17 1471.12 21851.87 -3473.23 16789.09 -752.03 1.85 -2.86 2.51 -3.28 -1.54 0.35 -1.93 -0.14 171.49 29.11 4.67 61.51
TC 999.35 2.68 867.57 129.09 0.40 -18.56 -30.48 -5.80 -15.25 -94.91 -301.21 115.29 -752.03 183.60 0.06 0.02 -0.16 0.32 0.14 0.00 0.17 0.02 -8.94 -1.05 1.79 -4.55
totFlat -0.79 -0.02 -0.87 0.09 0.00 0.03 0.04 0.03 0.02 -0.64 5.33 -1.76 1.85 0.06 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.01 -0.01 0.01
attFlat 1.26 0.02 1.19 0.05 0.00 0.01 -0.03 0.02 0.02 0.24 -5.54 1.22 -2.86 0.02 0.00 0.01 0.00 0.00 0.00 0.00 0.00 0.00 -0.05 -0.01 0.00 -0.01
decFlat -2.41 -0.02 -2.28 -0.10 -0.01 0.07 0.10 0.04 0.06 -0.75 6.07 -2.51 2.51 -0.16 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.06 0.01 -0.01 0.02
earlyDecayFlat 2.57 0.01 2.27 0.28 0.00 -0.02 -0.06 0.01 -0.01 -0.58 -4.17 0.47 -3.28 0.32 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.05 -0.01 0.00 -0.02
totSpecFlat 1.05 0.01 0.94 0.11 0.00 -0.03 -0.05 -0.02 -0.02 0.58 -1.74 1.42 -1.54 0.14 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.02 0.00 0.00 -0.01
attSpecFlat -0.17 0.00 -0.17 0.01 0.00 0.00 0.00 0.00 -0.01 0.17 1.92 0.13 0.35 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.00 0.00 0.00
decSpecFlat 1.31 0.01 1.16 0.13 0.00 -0.04 -0.06 -0.03 -0.03 0.72 -2.23 1.80 -1.93 0.17 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.02 0.00 0.00 -0.01
earlyDecaySpecFlat 0.12 0.00 0.11 0.02 0.00 -0.01 -0.01 0.00 -0.01 0.12 -0.10 0.24 -0.14 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
totCrest -90.78 -1.26 -82.96 -6.57 -0.26 -0.25 2.08 -1.13 -0.64 9.18 314.17 -53.07 171.49 -8.94 0.04 -0.05 0.06 -0.05 -0.02 0.01 -0.02 0.00 2.99 0.50 -0.04 0.88
attCrest -13.10 -0.26 -12.27 -0.57 -0.05 0.01 0.38 -0.10 -0.06 -3.27 49.56 -16.77 29.11 -1.05 0.01 -0.01 0.01 -0.01 0.00 0.00 0.00 0.00 0.50 0.11 -0.01 0.14
decCrest 27.04 -0.03 25.25 1.83 -0.01 -0.77 -0.61 -0.53 -0.63 1.29 -7.49 -0.43 4.67 1.79 -0.01 0.00 -0.01 0.00 0.00 0.00 0.00 0.00 -0.04 -0.01 0.29 -0.04
earlyDecayCrest -39.39 -0.40 -35.30 -3.68 -0.08 0.40 1.09 0.06 0.24 5.57 83.34 -13.80 61.51 -4.55 0.01 -0.01 0.02 -0.02 -0.01 0.00 -0.01 0.00 0.88 0.14 -0.04 0.34
# heatmap(cov_mat)
# image(cov_mat)
cov_mat.z <- scale(cov_mat, center = T, scale = T)
cov_mat.z_melt <- reshape2::melt(cov_mat.z)
ggplot(cov_mat.z_melt, aes(Var1, Var2)) +
  geom_tile(aes(fill = value)) +
  scale_fill_gradient2(low = "red", mid = "orange", high = "yellow") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.title = element_blank())

2 Modeling

Sofia, on 2020-07-07: “The main hypothesis is that playing instruction (N/C) will affect the stroke in a way that is perceivable”.

2.1 Univariate Models – transSC

Looking at the variables agreed upon I have decided to use as response variables in the regression models—implying that this is by no means set in stone—, it seems like a skewed normal link function would be appropriate to model them.

Let’s start with transSCtotDur. Using an extended model description language (Bates 2010; Bürkner 2018), going back to Wilkinson and Rogers’s (1973) modeling language, we write:

(m0_form <-bf(transSC ~ 1 + (1 | subj)))
transSC ~ 1 + (1 | subj) 

which claims that transSC is explained by (‘~’) an intercept, denoted by ‘1’, and an additional term ‘(1 | subj)’. The ‘1’ in parentheses again stands for the intercept, but the pipe ‘|’ assigns an intercept to each level of the factor ‘subj’. In this particular case this means that the model will estimate an individual intercept for each unique drummer listed in the data set column subj. These individual, or varying, intercepts are then used in informing the estimation of the population intercept.

Note: set MODEL to TRUE at the top of the script if you haven’t compiled/built your model yet.

if (MODEL) {
  m0 <- brm(m0_form,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m0 <- add_criterion(m0, 
                      "loo",
                      reloo = TRUE)
  save(m0, 
       file = "m0.rda")
} else {
  load("m0.rda")
}

This null model is also termed unconditional model because it has no grouping structure apart from individuals–the ‘subj’ bit in the model equation above. There is some variation in every natural data set. To make sure, it’s not just variation caused by different participants, we can calculate the intra-class correlation coefficient (ICC).

m0_icc <- ICC(m0, "subj")

The Null model’s ICC amounts to 0.76, which suggests that approx. 76 percent of the variation in the data set can be attributed to (or explained by) the grouping structure. This is highly unfortunate, as it does not leave a lot of variation to be explained by independent factors like instruction, or arm. The most likely reason for this high ICC value is the small sample size combined with high inter-individual variation. Small sample sizes combined with large trial numbers are less of a problem when subjects respond, on average, close to the population mean, even with large spread due to fluctuating alertness, increasing fatigue etc., .e.g. in reaction time paradigms. But here, with large intra- and inter-individual variation, this might become a problem.


Tab. 1. Model summary.

(m0_summary <- summary(m0, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ 1 + (1 | subj) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ student_t(3, 0, 136.5)
<lower=0> sigma ~ student_t(3, 0, 136.5)

Multilevel Hyperparameters:
~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)   116.78     34.58    69.10   201.98 1.01      690      699

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   747.16     39.66   665.75   828.48 1.01      703      962

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    64.98      1.50    62.02    68.02 1.00     2270     2035
alpha     3.14      0.33     2.53     3.85 1.00     1849     1521

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept in Tab. 1 amounts to roughly 747. In this model specification, the intercept is identical to the data set average. The table also shows that the inter-individual standard deviation (sd(Intercept)) is large compared to the unexplained variation (sigma). This led to the large ICC value above. By adding more and more independent factors to the model specification, we will later try to decrease \(\sigma\), i.e. ‘explain away’ as much remaining variation as possible.

pp_check(m0, ndraws = 100)

Fig. 12. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


The posterior predictive plot in Fig. 12 gives an impression on how the null model would generate data, given the parameters it estimated from the empirical data. As the thick line deviates from the modeled thin lines, especially on the left and the right side of the peak of the distribution, it is apparent that the model can be improved.

# conditional_effects(m0) 

2.1.1 Group Model

In this model, the intercept is complemented with second ‘main’ or population effect, the grouping variable patient:

(m1_form <- bf(transSC ~ 1 + patient + 
                 (1 | subj)))
transSC ~ 1 + patient + (1 | subj) 

I could have also specified the model without the explicit ‘1 +’, as the intercept is implicitly included in the model unless I explicitly exclude it. From now on I will always save some extra typing by refraining from explicitly indluding the intercept in models.

So now the model not only contains the individuals as grouping structure to ‘explain away’ variation, but also whether they belong to the patients or the healthy subjects.

The prior distributions in Bayesian models reflect the knowledge about the estimated parameters. The package brms automatically places weakly informative priors on parameters as soft contraints. But with more complex models involving varying parameter estimates, the statistical back end which does the heavy lifting, needs stronger priors particularly on these parameters, otherwise models don’t converge. The parameters most vulnerable to outliers in the data are the varying effects parameters, or random effects in traditional statistics. Therefore we place a stronger prior probability distribution over the estimate of the SD of individual intercepts:

(m1_prior <- set_prior("normal(0, 10)", class = "sd"))
sd ~ normal(0, 10)

and leave the rest of the priors as suggested by the package brms (see Tab. 2 for their priors).

if (MODEL) {
  m1 <- brm(m1_form,
            prior = m1_prior,
            init = "0",
            family = skew_normal(),
            data = drum_beats)
  m1 <- add_criterion(m1, 
                      "loo",
                      reloo = TRUE)
  save(m1, 
       file = "m1.rda")
} else {
  load("m1.rda")
}

Tab. 2. Model summary.

(m1_summary <- summary(m1, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ 1 + patient + (1 | subj) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 4)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 10)
<lower=0> sigma ~ student_t(3, 0, 136.5)

Multilevel Hyperparameters:
~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    49.11      4.74    40.50    59.08 1.00     2105     2383

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   762.85     24.82   715.40   812.67 1.00     1340     1601
patient1    -26.17     34.41   -93.26    40.92 1.00     1352     1536

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    65.12      1.53    62.26    68.27 1.00     3256     2760
alpha     3.26      0.33     2.64     3.94 1.00     3283     2630

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept in Tab. 2 amounts to roughly 763. In this model specification, the intercept is identical to the controlled strokes, while the patient1 value is the mean of the posterior distribution difference between the healthy subjects and the patients (-26). The table also shows that the interindividual standard deviation (sd(Intercept) \(\approx\) 49 is still large compared to the unexplained variation (sigma \(\approx\) 65).

pp_check(m1, ndraws = 100)

Fig. 13. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


The posterior predictive plot in Fig. 13 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the groups.

conditional_effects(m1,
                    dpar = "mu")

2.1.2 Group + Condition Model

In this model, the intercept is complemented with second ‘main’ or population effect, the instruction condition cond:

(m2_form <-bf(transSC ~ patient + cond + 
                (1 | subj)))
transSC ~ patient + cond + (1 | subj) 

So now the model also contains the individuals as grouping structure to ‘explain away’ variation, but also the manipulations.

The prior distributions in Bayesian models reflect the knowledge about the estimated parameters. The package brms automatically places weakly informative priors on parameters as soft contraints. But with more complex models involving varying parameter estimates, the statistical back end which does the heavy lifting, needs stronger priors particularly on these parameters, otherwise models don’t converge. The parameters most vulnerable to outliers in the data are the varying effects parameters, or random effects in traditional statistics. Therefore we place a stronger prior probability distribution over the estimate of the SD of individual intercepts:

(m2_prior <- c(set_prior("normal(0, 3)",   class = "sd"),
               set_prior("normal(0, 3)",   class = "sigma"),
               set_prior("normal(0, 2)",   class = "alpha")))
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user

and leave the rest of the priors as suggested by the package brms (see Tab. 3 for their priors).

if (MODEL) {
  m2 <- brm(m2_form,
            prior = m2_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m2 <- add_criterion(m2, 
                      "loo",
                      reloo = TRUE)
  save(m2, 
       file = "m2.rda")
} else {
  load("m2.rda")
}

Tab. 3. Model summary.

(m2_summary <- summary(m2, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient + cond + (1 | subj) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.25      1.51    26.45    32.24 1.00     2829     2544

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   779.96     15.20   750.21   808.91 1.01     1477     2027
patient1    -27.00     20.87   -67.13    14.96 1.01     1390     2000
condN       -41.31      3.43   -47.98   -34.63 1.00     4054     3059

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.37      0.98    50.48    54.27 1.00     4593     2495
alpha     1.80      0.24     1.36     2.28 1.00     4153     2874

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept in Tab. 3 amounts to roughly 780. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-41). The table also shows that the interindividual standard deviation (sd(Intercept) \(\approx\) 29 is still large compared to the unexplained variation (sigma \(\approx\) 52).

pp_check(m2, ndraws = 100)

Fig. 14. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


The posterior predictive plot in Fig. 14 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.

beautify_my_plot(plot(conditional_effects(m2,
                    dpar = "mu"),
                    plot = FALSE))

2.1.3 Group + ConditionVI Model

The last model included the instruction condition, more realistically reflecting the true structure of the data set. But it did not acknowledge that each subject executed several strokes in each of these conditions. This model includes a term with a varying intercept (VI) for condition to reflect just that, which will also assist in more realistically estimate the population effect of cond:

(m3_form <- bf(transSC ~ patient + cond + 
                 (1 | subj) + 
                 (1 | cond)))
transSC ~ patient + cond + (1 | subj) + (1 | cond) 

Including more varying parameters in the model requires the prior on them to be even stronger:

(m3_prior <- m2_prior)
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user
if (MODEL) {
  m3 <- brm(m3_form,
            prior = m3_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m3 <- add_criterion(m3, 
                      "loo",
                      reloo = TRUE)
  save(m3, 
       file = "m3.rda")
} else {
  load("m3.rda")
}

Tab. 4. Model summary.

(m3_summary <- summary(m3, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient + cond + (1 | subj) + (1 | cond) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.42      1.79     0.10     6.62 1.00     2230     1513

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.27      1.50    26.47    32.18 1.00     2628     2689

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   779.20     14.91   750.33   808.74 1.00     1685     2108
patient1    -25.86     20.09   -65.66    13.54 1.00     1704     2393
condN       -41.24      5.40   -52.30   -30.02 1.00     2838     2301

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.36      0.93    50.65    54.26 1.00     4618     2948
alpha     1.80      0.23     1.35     2.26 1.00     5002     2914

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept in Tab. 4 amounts to roughly 779. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-41). The table also shows that the inter-individual standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared to the unexplained variation (sigma \(\approx\) 52).

pp_check(m3, ndraws = 100)

Fig. 15. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


The posterior predictive plot in Fig. 15 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.

beautify_my_plot(plot(conditional_effects(m3,
                    dpar = "mu"),
                    plot = FALSE))

2.1.4 Group x ConditionVI

The last model included the instruction condition, more realistically reflecting the true structure of the data set. But it did not acknowledge that each subject executed several strokes in each of these conditions. This model includes a term with a varying intercept (VI) for condition to reflect just that, which will also assist in more realistically estimate the population effect of cond:

(m4_form <- bf(transSC ~ patient * cond + 
                 (1 | subj) + 
                 (1 | cond)))
transSC ~ patient * cond + (1 | subj) + (1 | cond) 

Including more varying parameters in the model requires the prior on them to be even stronger:

(m4_prior <- m3_prior)
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user
if (MODEL) {
  m4 <- brm(m4_form,
            prior = m4_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m4 <- add_criterion(m4, 
                      "loo",
                      reloo = TRUE)
  save(m4, 
       file = "m4.rda")
} else {
  load("m4.rda")
}

Tab. 5. Model summary.

(m4_summary <- summary(m4, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient * cond + (1 | subj) + (1 | cond) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.43      1.84     0.08     6.90 1.00     2489     1546

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.23      1.48    26.36    32.17 1.00     3437     3199

Regression Coefficients:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        780.51     15.43   751.24   811.37 1.00     1548     2068
patient1         -28.11     21.00   -70.61    11.97 1.00     1554     2272
condN            -42.62      6.34   -55.38   -29.58 1.00     2787     2530
patient1:condN     2.44      6.33    -9.93    14.50 1.00     5037     3014

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.39      0.91    50.66    54.21 1.00     5152     3046
alpha     1.81      0.24     1.36     2.31 1.00     4824     3033

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

The intercept in Tab. 5 amounts to roughly 781. In this model specification, the intercept is identical to the controlled strokes, while the condN value is the mean of the posterior distribution difference between controlled and normal strokes (-43). The table also shows that the inter-individual standard deviation (sd(Intercept) \(\approx\) 29 is not large anymore compared to the unexplained variation (sigma \(\approx\) 52).

pp_check(m4, ndraws = 100)

Fig. 16. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


The posterior predictive plot in Fig. 16 gives an impression on how this model would generate data, given the parameters it estimated from the empirical data. There is not much change from Fig. 12, which is not very surprising given the small estimated difference between the conditions.

beautify_my_plot(plot(conditional_effects(m4,
                    dpar = "mu"),
                    plot = FALSE))

2.1.5 Group x ConditionVI + Arm Model

(m5_form <- bf(transSC ~ patient * cond + arm +
                 (1 | subj) + 
                 (1 | cond)))
transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) 
(m5_prior <- m4_prior)
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user
if (MODEL) {
  m5 <- brm(m5_form,
            prior = m5_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m5 <- add_criterion(m5, 
                      "loo",
                      reloo = TRUE)
  save(m5, 
       file = "m5.rda")
} else {
  load("m5.rda")
}

Tab. 6. Model summary.

(m5_summary <- summary(m5, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.38      1.81     0.10     6.85 1.00     2808     1452

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.26      1.48    26.54    32.26 1.00     3002     2712

Regression Coefficients:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        778.69     15.22   750.01   808.84 1.00     1599     2253
patient1         -27.48     21.16   -67.63    13.57 1.00     1482     1997
condN            -42.50      6.19   -54.84   -30.18 1.00     2940     2333
armND              1.94      3.17    -4.19     8.16 1.00     4885     2750
patient1:condN     2.44      6.31    -9.70    15.13 1.00     3896     2989

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.39      0.93    50.62    54.24 1.00     4932     3065
alpha     1.80      0.23     1.36     2.26 1.00     4264     2854

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Population-Level Effects section of Tab. 6 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.

pp_check(m5, ndraws = 100)

Fig. 17. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m5,
                    dpar = "mu"),
                    plot = FALSE))

2.1.6 Group x Condition + ArmVI Model

(m6_form <-bf(transSC ~ patient * cond + arm +
                (1 | subj) + 
                (1 | cond) + 
                (1 | arm)))
transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm) 
(m6_prior <- m5_prior)
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user
if (MODEL) {
  m6 <- brm(m6_form,
            prior = m6_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m6 <- add_criterion(m6, 
                      "loo",
                      reloo = TRUE)
  save(m6, 
       file = "m6.rda")
} else {
  load("m6.rda")
}

Tab. 7. Model summary.

(m6_summary <- summary(m6, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient * cond + arm + (1 | subj) + (1 | cond) + (1 | arm) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.36      1.80     0.09     6.72 1.00     1732     1395

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.37      1.78     0.09     6.61 1.00     2597     1930

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.27      1.51    26.41    32.38 1.00     3277     2942

Regression Coefficients:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept        778.89     16.31   746.65   810.19 1.00     1636     1936
patient1         -27.71     21.59   -69.78    14.20 1.00     1737     1690
condN            -42.73      6.05   -55.41   -30.77 1.00     2625     2270
armND              1.96      5.20    -8.77    12.69 1.00     2936     2231
patient1:condN     2.57      6.17    -9.32    14.65 1.00     4105     2772

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.36      0.93    50.50    54.21 1.00     4749     2764
alpha     1.80      0.23     1.35     2.27 1.00     5175     3159

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Population-Level Effects section of Tab. 7 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.

pp_check(m6, ndraws = 100)

Fig. 18. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m6,
                    dpar = "mu"),
                    plot = FALSE))

2.1.7 Group x Condition x ArmVI Model

(m6a_form <-bf(transSC ~ patient * cond * arm +
                 (1 | subj) + 
                 (1 | cond) + 
                 (1 | arm)))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
(m6a_prior <- m6_prior)
        prior class coef group resp dpar nlpar   lb   ub source
 normal(0, 3)    sd                            <NA> <NA>   user
 normal(0, 3) sigma                            <NA> <NA>   user
 normal(0, 2) alpha                            <NA> <NA>   user
if (MODEL) {
  m6a <- brm(m6a_form,
             prior = m6a_prior,
             family = skew_normal(),
             init = "0",
             data = drum_beats)
  m6a <- add_criterion(m6a, 
                       "loo",
                       reloo = TRUE)
  save(m6a, 
       file = "m6a.rda")
} else {
  load("m6a.rda")
}

Tab. 8. Model summary.

(m6a_summary <- summary(m6a, 
                        priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
Intercept ~ student_t(3, 735, 136.5)
<lower=0> sd ~ normal(0, 3)
<lower=0> sigma ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.38      1.82     0.07     6.83 1.00     2462     1425

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.44      1.84     0.10     6.72 1.00     2523     1781

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    29.39      1.47    26.64    32.40 1.00     3165     2828

Regression Coefficients:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept              784.68     15.62   754.83   815.57 1.00     1612
patient1               -42.15     21.48   -83.41    -0.29 1.00     1582
condN                  -42.56      7.49   -57.33   -27.57 1.00     2083
armND                   -7.65      7.22   -22.05     6.73 1.00     2220
patient1:condN           9.10      8.96    -8.29    26.78 1.00     2454
patient1:armND          26.71      8.44    10.13    43.50 1.00     2613
condN:armND             -0.19      8.56   -17.05    16.61 1.00     2081
patient1:condN:armND   -13.95     11.97   -38.00     9.13 1.00     2073
                     Tail_ESS
Intercept                2287
patient1                 2147
condN                    2445
armND                    2817
patient1:condN           2823
patient1:armND           2724
condN:armND              2508
patient1:condN:armND     2766

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    52.16      0.91    50.44    54.03 1.00     4903     3201
alpha     1.68      0.24     1.21     2.18 1.00     4802     2979

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Population-Level Effects section of Tab. 8 are now two estimates for the skewness parameter alpha (one for the controlled [alpha_Intercept] and one for the normal condition [alpha_condN]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on instruction condition.

pp_check(m6a, ndraws = 100)

Fig. 19. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m6a,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m6a,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                    plot = FALSE))

2.1.8 Modeling Condition-specific Variation

As is apparent from Fig. 1, the controlled condition yielded broader density distributions in most descriptors. Hence, we model condition-dependent variation in the next model:

(m7_form <-bf(transSC ~ patient * cond * arm + 
                (1 | subj) + 
                (1 | cond) +
                (1 | arm),
              sigma ~ cond))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ cond

This model has now two outcomes, not just one, as the models before. Spread is estimated as sigma, and it varies conditional on instruction condition.

Here’s the non-standard (additional) prior:

(m7_prior <- c(set_prior("normal(0, 5)", class = "sd"),
               set_prior("normal(0, 2)",  class = "alpha"),
               set_prior("normal(0, 10)", class = "Intercept", dpar = "sigma"),
               set_prior("normal(0, 10)", class = "b",         dpar = "sigma")
))
         prior     class coef group resp  dpar nlpar   lb   ub source
  normal(0, 5)        sd                             <NA> <NA>   user
  normal(0, 2)     alpha                             <NA> <NA>   user
 normal(0, 10) Intercept                 sigma       <NA> <NA>   user
 normal(0, 10)         b                 sigma       <NA> <NA>   user
if (MODEL) {
  m7 <- brm(m7_form,
            prior = m7_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m7 <- add_criterion(m7, 
                      "loo",
                      reloo = TRUE)
  save(m7, 
       file = "m7.rda")
} else {
  load("m7.rda")
}

Tab. 9. Model summary.

(m7_summary <- summary(m7, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ cond
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 5)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     4.01      3.05     0.13    11.09 1.00     2457     1648

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     3.97      2.97     0.15    11.22 1.00     3501     2247

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    35.46      2.38    30.98    40.35 1.00     4063     3041

Regression Coefficients:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept              788.91     20.56   747.69   830.07 1.00     1641
sigma_Intercept          4.30      0.04     4.23     4.38 1.00     5439
patient1               -41.29     26.62   -93.30    11.66 1.00     1621
condN                  -45.64      9.78   -64.94   -26.17 1.00     2365
armND                   -9.12     10.98   -30.01    13.56 1.00     2578
patient1:condN           7.46     10.14   -11.97    27.61 1.00     3092
patient1:armND          24.81     11.66     2.43    48.21 1.00     2642
condN:armND              0.11      9.70   -18.49    18.96 1.00     2610
patient1:condN:armND   -10.28     13.92   -37.46    17.60 1.00     2386
sigma_condN             -0.46      0.06    -0.57    -0.34 1.00     4995
                     Tail_ESS
Intercept                2320
sigma_Intercept          2463
patient1                 2120
condN                    2489
armND                    1771
patient1:condN           2831
patient1:armND           2941
condN:armND              2884
patient1:condN:armND     2882
sigma_condN              2696

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
alpha     2.25      0.29     1.72     2.86 1.00     4782     3109

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Links section of the summary (Tab. 9) we note that sigma is no longer modeled on the identity scale but on the \(\log_{2}\) scale. In the Population-Level Effects section of Tab. 9 we now find not only the estimates for the Intercept and condN but also two estimates for sigma (one for the controlled [sigma_Intercept] and one for the normal condition [sigma_condN]). Because of this, there is no entry for sigma in the Family Specific section anymore. All this due to our explicit modeling of sigma conditional on instruction condition.

pp_check(m7, ndraws = 100)

Fig. 20. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m7,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m7,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m7,
                    dpar = "sigma"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m7,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                    plot = FALSE))

Fig. 21. Conditional plot.


In Fig. 21 we see a difference between the mean (‘mu’) estimated value for transSC depending on instruction condition (A), but no difference between dominant and non-dominant arm (B). Consequently, there is no interaction in (C).

The spread (sigma) of the estimated distributions also differs between conditions (D), but again not between arms (E). The latter is not surprising because we modeled sigma to vary conditional on condition, not arm.

2.1.9 Modeling Arm-specific Skewness

At least in some descriptors in Fig. 1 the skewness of the distribution seems to change depending on the arm. The skew normal distribution is a generalization of the Gaussian distribution, allowing for the additional shape-parameter skewness (asymmetry) to vary. Hence, we model this side-dependent skewness in the next model:

(m8_form <-bf(transSC ~ patient * cond * arm + 
                (1 | subj) + 
                (1 | cond) +
                (1 | arm),
              sigma ~ cond,
              alpha ~ arm))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ cond
alpha ~ arm
(m8_prior <- c(set_prior("normal(0, 3)",  class = "sd"),
               set_prior("normal(0, 10)", class = "Intercept", dpar = "sigma"),
               set_prior("normal(0, 10)", class = "b",         dpar = "sigma"),
               set_prior("normal(0, 2)",  class = "Intercept", dpar = "alpha"),
               set_prior("normal(0, 2)",  class = "b",         dpar = "alpha"))
)
         prior     class coef group resp  dpar nlpar   lb   ub source
  normal(0, 3)        sd                             <NA> <NA>   user
 normal(0, 10) Intercept                 sigma       <NA> <NA>   user
 normal(0, 10)         b                 sigma       <NA> <NA>   user
  normal(0, 2) Intercept                 alpha       <NA> <NA>   user
  normal(0, 2)         b                 alpha       <NA> <NA>   user
if (MODEL) {
  m8 <- brm(m8_form,
            prior = m8_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m8 <- add_criterion(m8, 
                      "loo",
                      reloo = TRUE)
  save(m8, 
       file = "m8.rda")
} else {
  load("m8.rda")
}

Tab. 10. Model summary.

(m8_summary <- summary(m8, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ cond
         alpha ~ arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.35      1.78     0.09     6.59 1.00     2975     1680

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.42      1.83     0.10     6.75 1.00     3645     2387

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    27.38      1.51    24.49    30.49 1.00     2970     2278

Regression Coefficients:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept              787.52     15.14   757.90   816.85 1.00     1551
sigma_Intercept          4.33      0.04     4.26     4.41 1.00     4488
alpha_Intercept          1.69      0.41     0.84     2.40 1.00     4289
patient1               -43.17     20.57   -82.05    -1.86 1.00     1674
condN                  -44.43      8.36   -61.25   -28.21 1.00     2119
armND                   -9.63      9.07   -26.73     8.35 1.00     2117
patient1:condN           9.38     10.65   -11.30    30.57 1.00     2068
patient1:armND          27.49     11.77     3.94    50.13 1.00     1892
condN:armND             -1.48      9.57   -20.15    17.16 1.00     1962
patient1:condN:armND    -9.96     13.91   -37.79    17.02 1.00     1881
sigma_condN             -0.51      0.06    -0.62    -0.40 1.00     4226
alpha_armND              1.68      0.54     0.67     2.80 1.00     3751
                     Tail_ESS
Intercept                2375
sigma_Intercept          2913
alpha_Intercept          2079
patient1                 2309
condN                    2255
armND                    2509
patient1:condN           2545
patient1:armND           2594
condN:armND              2428
patient1:condN:armND     2228
sigma_condN              3109
alpha_armND              2446

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Population-Level Effects section of Tab. 10 are now two estimates for the skewness parameter alpha (one for the dominant [alpha_Intercept] and one for the non-dominant arm [alpha_armND]), while there is no entry for alpha in the Family Specific section anymore. All this due to our explicit modeling of alpha conditional on arm.

pp_check(m8, ndraws = 100)

Fig. 22. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "mu"),
                      plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "sigma"),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "alpha"),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m8,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 23. Conditional plots.


In Fig. 23 there is a clear difference in estimated mean (‘mu’) of \(transSC\) conditional on instruction (A), but not on arm (B). Consequently, there is no interaction between instruction and arm (C). On the other hand, there is a difference in estimated skewness (‘alpha’) between sides (E), but not instructions (D). (F) follows from that. ‘sigma’ shows a clear difference between instructions (G), not so for side (H), which is mirrored in (I).

2.1.10 Modeling Influence of Arm on Shape

To be able to see whether skewness and spread are depending on both experimental manipulation and arm, the model formulas for ‘sigma’ and ‘alpha’ are updated accordingly:

(m9_form <-bf(transSC ~ patient * cond * arm + 
                (1 | subj) + 
                (1 | cond) +
                (1 | arm),
              sigma ~ cond,
              alpha ~ arm))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ cond
alpha ~ arm
(m9_prior <- m8_prior)
         prior     class coef group resp  dpar nlpar   lb   ub source
  normal(0, 3)        sd                             <NA> <NA>   user
 normal(0, 10) Intercept                 sigma       <NA> <NA>   user
 normal(0, 10)         b                 sigma       <NA> <NA>   user
  normal(0, 2) Intercept                 alpha       <NA> <NA>   user
  normal(0, 2)         b                 alpha       <NA> <NA>   user
if (MODEL) {
  m9 <- brm(m9_form,
            prior = m9_prior,
            family = skew_normal(),
            init = "0",
            data = drum_beats)
  m9 <- add_criterion(m9, 
                      "loo",
                      reloo = TRUE)
  save(m9, 
       file = "m9.rda")
} else {
  load("m9.rda")
}

Tab. 11. Model summary.

(m9_summary <- summary(m9, 
                       priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ cond
         alpha ~ arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.41      1.80     0.10     6.62 1.00     3102     1803

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.40      1.81     0.12     6.79 1.00     3362     2082

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    27.40      1.51    24.51    30.42 1.00     3438     3097

Regression Coefficients:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept              787.88     15.43   757.82   819.11 1.00     1707
sigma_Intercept          4.33      0.04     4.26     4.41 1.00     4444
alpha_Intercept          1.70      0.40     0.89     2.45 1.00     4143
patient1               -44.33     21.11   -86.36    -2.99 1.00     1676
condN                  -44.87      8.20   -61.12   -28.53 1.00     2114
armND                  -10.20      9.12   -27.85     8.26 1.00     1934
patient1:condN          10.29     10.46   -10.00    31.03 1.00     1816
patient1:armND          28.35     11.81     5.73    51.76 1.00     1713
condN:armND             -0.81      9.48   -19.85    17.70 1.00     1904
patient1:condN:armND   -11.21     13.77   -38.63    15.43 1.00     1541
sigma_condN             -0.52      0.06    -0.63    -0.40 1.00     4153
alpha_armND              1.66      0.55     0.60     2.81 1.00     3858
                     Tail_ESS
Intercept                2215
sigma_Intercept          2965
alpha_Intercept          2447
patient1                 2169
condN                    2946
armND                    2731
patient1:condN           2771
patient1:armND           2531
condN:armND              2694
patient1:condN:armND     2299
sigma_condN              3089
alpha_armND              2679

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

In the Population-Level Effects section of Tab. 11 are now three estimates for the skewness parameter alpha, as well as three estimates for spread (‘sigma’).

pp_check(m9, ndraws = 100)

Fig. 24. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m9,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m9,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m9,
                    dpar = "sigma"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m9,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m9,
                    dpar = "alpha"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m9,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 25. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.


In Fig. 25 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.

There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).

The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).

There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).

From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.

2.1.11 Modeling Interactions of Arm and Instruction on Shape and Spread

If an interaction between cond and arm should be part of the model explaining the outcome transSC, then this interactin should be also present in the part of the model explaining sigma and alpha:

(m10_form <-bf(transSC ~ patient * cond * arm + 
                 (1 | subj) + 
                 (1 | cond) +
                 (1 | arm),
               sigma ~ cond * arm,
               alpha ~ cond * arm))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ cond * arm
alpha ~ cond * arm
(m10_prior <- m9_prior)
         prior     class coef group resp  dpar nlpar   lb   ub source
  normal(0, 3)        sd                             <NA> <NA>   user
 normal(0, 10) Intercept                 sigma       <NA> <NA>   user
 normal(0, 10)         b                 sigma       <NA> <NA>   user
  normal(0, 2) Intercept                 alpha       <NA> <NA>   user
  normal(0, 2)         b                 alpha       <NA> <NA>   user
if (MODEL) {
  m10 <- brm(m10_form,
             prior = m10_prior,
             family = skew_normal(),
             init = "0",
             data = drum_beats)
  m10 <- add_criterion(m10, 
                       "loo",
                       reloo = TRUE)
  save(m10, 
       file = "m10.rda")
} else {
  load("m10.rda")
}

Tab. 12. Model summary.

(m10_summary <- summary(m10, 
                        priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ cond * arm
         alpha ~ cond * arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.41      1.83     0.13     6.86 1.00     3683     2088

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.44      1.79     0.10     6.69 1.00     2636     1699

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    27.40      1.50    24.65    30.46 1.00     3648     2642

Regression Coefficients:
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept              787.53     15.38   757.35   817.04 1.00     1984
sigma_Intercept          4.29      0.05     4.20     4.40 1.00     3976
alpha_Intercept          4.40      0.65     3.23     5.72 1.00     4080
patient1               -38.10     20.72   -78.43     2.13 1.00     2100
condN                  -44.27      7.81   -59.75   -28.73 1.00     2882
armND                   -9.11      9.43   -27.43     9.42 1.00     2618
patient1:condN           2.15      9.19   -16.33    19.86 1.00     2719
patient1:armND          31.16     10.83     9.38    52.14 1.00     2578
condN:armND             -2.98      9.90   -22.31    16.33 1.00     2465
patient1:condN:armND   -14.15     13.04   -39.40    11.65 1.00     2271
sigma_condN             -0.52      0.07    -0.66    -0.38 1.00     3888
sigma_armND              0.13      0.07    -0.00     0.27 1.00     3423
sigma_condN:armND       -0.18      0.09    -0.36    -0.00 1.00     3469
alpha_condN             -5.35      0.87    -6.92    -3.46 1.00     2477
alpha_armND              2.05      1.03     0.16     4.16 1.00     3322
alpha_condN:armND        1.05      1.16    -1.35     3.24 1.00     2460
                     Tail_ESS
Intercept                2510
sigma_Intercept          3305
alpha_Intercept          2860
patient1                 2414
condN                    2815
armND                    3087
patient1:condN           2812
patient1:armND           2815
condN:armND              2630
patient1:condN:armND     2618
sigma_condN              3186
sigma_armND              2930
sigma_condN:armND        2914
alpha_condN              2856
alpha_armND              2450
alpha_condN:armND        2580

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(m10, ndraws = 100)

Fig. 26. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m10,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m10,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m10,
                    dpar = "sigma"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m10,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m10,
                    dpar = "alpha"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m10,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 27. Conditional plots. Had to take out all the fancy figure polishing. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise. Instead, they are now numbered from A to F, plus one un-numbered plot with the three-way interaction, for each of the estimated outcomes transSC, sigma, and alpha


In Fig. 27 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.

There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).

The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).

There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).

From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.

2.1.12 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread

(m11_form <-bf(transSC ~ patient * cond * arm + 
                 (1 | subj) + 
                 (1 | cond) +
                 (1 | arm),
               sigma ~ patient * cond * arm,
               alpha ~ patient * cond * arm))
transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
(m11_prior <- m10_prior)
         prior     class coef group resp  dpar nlpar   lb   ub source
  normal(0, 3)        sd                             <NA> <NA>   user
 normal(0, 10) Intercept                 sigma       <NA> <NA>   user
 normal(0, 10)         b                 sigma       <NA> <NA>   user
  normal(0, 2) Intercept                 alpha       <NA> <NA>   user
  normal(0, 2)         b                 alpha       <NA> <NA>   user
if (MODEL) {
  m11 <- brm(m11_form,
             prior = m11_prior,
             family = skew_normal(),
             init = "0",
             data = drum_beats)
  m11 <- add_criterion(m11, 
                       "loo",
                       reloo = TRUE)
  save(m11, 
       file = "m11.rda")
} else {
  load("m11.rda")
}

Tab. 13. Model summary.

(m11_summary <- summary(m11, 
                        priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transSC ~ patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ patient * cond * arm
         alpha ~ patient * cond * arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 10)
Intercept ~ student_t(3, 735, 136.5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 10)
<lower=0> sd ~ normal(0, 3)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.46      1.85     0.12     6.94 1.00     4068     2366

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     2.40      1.81     0.12     6.90 1.00     4364     2706

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)    27.31      1.49    24.51    30.26 1.00     6072     3196

Regression Coefficients:
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                    790.58     15.68   759.87   820.68 1.00     1928
sigma_Intercept                4.34      0.07     4.21     4.47 1.00     2907
alpha_Intercept                4.42      0.77     3.00     6.00 1.00     4572
patient1                     -46.81     21.23   -88.55    -5.94 1.00     2013
condN                        -47.05      8.30   -63.05   -30.07 1.00     2747
armND                        -12.06      9.99   -31.54     7.96 1.00     2188
patient1:condN                10.77      9.98    -9.04    30.28 1.00     2459
patient1:armND                41.45     12.83    16.20    66.79 1.00     2070
condN:armND                    1.32     10.43   -19.04    21.39 1.00     2085
patient1:condN:armND         -26.89     14.85   -55.77     2.18 1.00     1957
sigma_patient1                -0.15      0.10    -0.35     0.05 1.00     2531
sigma_condN                   -0.54      0.11    -0.74    -0.32 1.00     3088
sigma_armND                    0.05      0.10    -0.16     0.23 1.00     2477
sigma_patient1:condN           0.11      0.15    -0.18     0.40 1.00     2678
sigma_patient1:armND           0.21      0.14    -0.06     0.49 1.00     2243
sigma_condN:armND             -0.03      0.13    -0.29     0.22 1.00     2595
sigma_patient1:condN:armND    -0.33      0.19    -0.69     0.04 1.00     2537
alpha_patient1                -0.61      1.03    -2.64     1.46 1.00     4839
alpha_condN                   -5.07      1.10    -7.15    -2.98 1.00     2807
alpha_armND                    0.64      1.20    -1.60     3.03 1.00     4209
alpha_patient1:condN           0.49      1.21    -1.86     2.88 1.00     3406
alpha_patient1:armND           2.95      1.32     0.41     5.53 1.00     4069
alpha_condN:armND              1.25      1.31    -1.35     3.80 1.00     3592
alpha_patient1:condN:armND    -1.36      1.37    -4.08     1.36 1.00     4208
                           Tail_ESS
Intercept                      2519
sigma_Intercept                3364
alpha_Intercept                3149
patient1                       2584
condN                          2852
armND                          2727
patient1:condN                 3169
patient1:armND                 2793
condN:armND                    2875
patient1:condN:armND           2639
sigma_patient1                 3017
sigma_condN                    3142
sigma_armND                    2916
sigma_patient1:condN           2823
sigma_patient1:armND           2913
sigma_condN:armND              3185
sigma_patient1:condN:armND     2737
alpha_patient1                 2941
alpha_condN                    3011
alpha_armND                    3205
alpha_patient1:condN           3008
alpha_patient1:armND           3298
alpha_condN:armND              3017
alpha_patient1:condN:armND     3177

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(m11, ndraws = 100)

Fig. 28. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m11,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m11,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11,
                    dpar = "sigma"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11,
                    dpar = "alpha"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 29. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.


In Fig. 29 the dots are the means of the posterior distributions of the respective estimate, and can be interpreted similarly to the point estimates in traditional frequentist statistics. The error bars represent 95% credible intervals, which are interpreted as comprising the true value with 95% probability, given the model and the data.

There is a difference in estimated the mean (‘mu’) of \(transSC\) conditional on group (A) and condition (B), not so on arm (C). Likewise, there is a two-way interaction between group and arm (D), but not between condition and arm (E). Consequently, there is no three-way interaction (F).

The remaining, unexplained variation ‘sigma’ shows a clear difference between instructions (H) and lesser so for side (I), but albeit their opposing trends, there is no interaction (J thru M).

There is a difference in estimated skewness (‘alpha’) both between instruction (O) and side (P), but with opposing trends. Nevertheless, there is no interaction between them (Q thru T).

From a data-driven perspective this should lead to a simplification of the model, leaving out all the unnecessary interaction terms. A theory-driven approach though requires the final model to be ‘full’ in the sense that all explanatory variables must interact with all others, implying that not finding an interaction between factors either is due to there really not being an interaction–rendering the model mis-specified and proving the underlying theory wrong–, or due to insufficient data. Since the data set at hand is very small I vote to blame the failure to find interactions between all factors on the small sample size and would therefore keep all interactions in the model.

2.1.13 Model Comparison

We compare models by their estimated log-posterior density (elpd). The smaller this value the better a model predicts the data, despite the penalty for additional covariates. When the difference between two models is more than two SE apart they are considered to be ‘different enough’ to warrant acceptance of one over the other despite possibly smaller parsimony.

univ_model_compar <- loo_compare(m0, m2, m3, m5, m6, m6a, m7, m8, m9, m10, m11)
print(univ_model_compar, simplify = T)
    elpd_diff se_diff
m10    0.0       0.0 
m11   -2.1       3.7 
m9   -40.7       9.0 
m8   -40.8       8.9 
m7   -43.6      10.0 
m6a  -83.1      15.8 
m3   -83.2      16.7 
m2   -83.5      16.7 
m6   -85.1      16.8 
m5   -85.4      16.8 
m0  -121.7      15.5 

The model with the smallest elpd is m10, which corresponds to the second to last model (see Modeling Interactions of Arm and Instruction on Shape and Spread). The difference between the best and the second-best model (m11, [Modeling the Influence of Group on Shape and Spread]) is -2.1, which is less than roughly two SE (2 x 3.67 = 7.34) away from 0. Therefore, we consider the best fitting and the second-best fitting models (roughly) equivalent. Based on theoretical grounds, modeling the three estimates mu, sigma and alpha should be done with the full model–i.e., containing all the main effect interactions–if there is no risk of overfitting. Therefore, in the following we use the full model ([Modeling the Influence of Group on Shape and Spread]).

2.1.14 Conclusion

The fact that the full model turned out to be the second-best suggests (see above) that the data really are best fit by this model, and that the patient term and the interaction terms really should remain in the model.

wrap_plots(plot(conditional_effects(m11,
                                    dpar = "mu",
                                    effects = "patient:cond",
                                    conditions = conditions),
                # facet_args = list(arm = "label_value"),
                plot = FALSE)[[1]] + 
             scale_color_colorblind() +
             theme(legend.position = "none"),
           plot(conditional_effects(m11,
                                    dpar = "sigma",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]] +
             scale_color_colorblind()  + 
             theme(legend.position = "none"), 
           plot(conditional_effects(m11,
                                    dpar = "alpha",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]]) +
  scale_color_colorblind() + 
  plot_annotation(tag_levels = "A")

Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the full model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).


Fig. 30A suggests that \(\mu\) is higher for healthy subjects than for patients, and higher for controlled strokes than for normal strokes which is really odd on an intuitive level and seems contradictive (at least to me). I would’ve expected transSC to be at some optimal level in healthy subjects while performing normal strokes with their dominant side, and decrease from there. Yet it does not seem to be so simple. The biggest and clearest difference between any of the experimental factors is between instruction conditions “normal” and “controlled” (see Tab. 13 for the credible intervals of these estimates). Of the interactions shown in Fig. 30A, only the estimate for patient1:armND is sufficiently far from zero, while the central 95% of the posterior distribution of patient1:condN:armND is, after all, located mostly to the left of zero.

Fig. 30B might clarify why this is the case, because it shows a higher unexplained variation in the data for controlled strokes as compared to normal strokes which was expected due to controlled strokes are executed with less automatized motor programs (at least that’s what the lay person at the keyboard thinks). And since higher variation is often accompanied by more extreme values, these extreme values might drive the mean of the posterior distribution of \(\mu\) up. Interestingly, \(\sigma\) is estimated to have almost the same low value across groups and sides in the normal instruction condition, whereas it is unanimously high in the controlled condition.

Fig. 30C gives the impression that in the normal condition the distribution of transSC is roughly symmetrical around zero, while in the controlled condition there is an increased amount of higher transSC values making the distribution heavy-tailed on the right side (see e.g., Wikipedia on skewness).

Summarizing these results, transSC is, on average, lower in normal strokes while approximately normally distributed, and with lower spread than in controlled strokes. The fact that there is a difference in the location parameter \(\mu\) between the two instructions suggests that this descriptor might serve as a cue in listening tests.

2.2 Univariate Models – transFlat

2.2.1 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread

(m11_transFlat_form <-bf(transFlat ~ 0 + Intercept + patient * cond * arm + 
                           (1 | subj) + 
                           (1 | cond) +
                           (1 | arm),
                         sigma ~ patient * cond * arm,
                         alpha ~ patient * cond * arm))
transFlat ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
(m11_transFlat_prior <- c(set_prior("normal(0.8, 10)", class = "b", coef = "Intercept"),
                          set_prior("normal(0, 10)",   class = "b"),
                          set_prior("normal(0, 5)",   class = "b", coef = "armND"),
                          set_prior("normal(0, 5)",   class = "b", coef = "condN"),
                          set_prior("normal(0, 0.01)", class = "sd"),
                          set_prior("normal(0, 5)",   class = "Intercept", dpar = "sigma"),
                          set_prior("normal(0, 5)",   class = "b",         dpar = "sigma"),
                          set_prior("normal(0, 2)",   class = "Intercept", dpar = "alpha"),
                          set_prior("normal(0, 2)",   class = "b",         dpar = "alpha")
))
           prior     class      coef group resp  dpar nlpar   lb   ub source
 normal(0.8, 10)         b Intercept                        <NA> <NA>   user
   normal(0, 10)         b                                  <NA> <NA>   user
    normal(0, 5)         b     armND                        <NA> <NA>   user
    normal(0, 5)         b     condN                        <NA> <NA>   user
 normal(0, 0.01)        sd                                  <NA> <NA>   user
    normal(0, 5) Intercept                      sigma       <NA> <NA>   user
    normal(0, 5)         b                      sigma       <NA> <NA>   user
    normal(0, 2) Intercept                      alpha       <NA> <NA>   user
    normal(0, 2)         b                      alpha       <NA> <NA>   user
if (MODEL) {
  m11_transFlat <- brm(m11_transFlat_form,
                       prior = m11_transFlat_prior,
                       family = skew_normal(),
                       init = "0",
                       data = drum_beats)
  m11_transFlat <- add_criterion(m11_transFlat, 
                                 "loo",
                                 reloo = TRUE)
  save(m11_transFlat, 
       file = "m11_transFlat.rda")
} else {
  load("m11_transFlat.rda")
}

Tab. 14. Model summary. The model has not converged.

(m11_transFlat_summary <- summary(m11_transFlat, 
                                  priors = TRUE))
Warning: There were 8 divergent transitions after warmup. Increasing
adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transFlat ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ patient * cond * arm
         alpha ~ patient * cond * arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b ~ normal(0, 10)
b_armND ~ normal(0, 5)
b_condN ~ normal(0, 5)
b_Intercept ~ normal(0.8, 10)
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 5)
<lower=0> sd ~ normal(0, 0.01)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.01      0.01     0.00     0.02 1.00     3240     2416

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.01      0.01     0.00     0.02 1.00     2715     1855

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.02      0.00     0.02     0.03 1.00     2888     2807

Regression Coefficients:
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept               -3.53      0.07    -3.65    -3.40 1.00     2679
alpha_Intercept               -4.26      0.74    -5.85    -2.93 1.00     5025
Intercept                      0.70      0.02     0.67     0.74 1.00     1817
patient1                       0.05      0.02     0.02     0.08 1.00     1548
condN                          0.02      0.01    -0.01     0.05 1.00     2225
armND                          0.01      0.01    -0.02     0.04 1.00     2364
patient1:condN                -0.01      0.00    -0.02    -0.00 1.00     3165
patient1:armND                -0.02      0.01    -0.03    -0.01 1.00     2333
condN:armND                    0.01      0.01     0.00     0.02 1.00     2554
patient1:condN:armND          -0.01      0.01    -0.02     0.00 1.00     2313
sigma_patient1                -0.22      0.10    -0.41    -0.04 1.00     2343
sigma_condN                   -0.15      0.10    -0.35     0.04 1.00     2417
sigma_armND                    0.45      0.09     0.28     0.63 1.00     2253
sigma_patient1:condN          -0.16      0.14    -0.44     0.11 1.00     2125
sigma_patient1:armND          -0.29      0.13    -0.53    -0.03 1.00     2128
sigma_condN:armND             -0.76      0.13    -1.01    -0.52 1.00     1857
sigma_patient1:condN:armND     0.48      0.18     0.12     0.84 1.00     1924
alpha_patient1                 0.03      1.05    -2.06     2.02 1.00     5137
alpha_condN                    2.37      0.91     0.65     4.18 1.00     4646
alpha_armND                   -1.38      1.04    -3.48     0.59 1.00     4454
alpha_patient1:condN           0.90      1.24    -1.55     3.34 1.00     4281
alpha_patient1:armND          -0.72      1.38    -3.47     1.97 1.00     5135
alpha_condN:armND              0.67      1.24    -1.71     3.12 1.00     4001
alpha_patient1:condN:armND     0.45      1.54    -2.51     3.50 1.00     5173
                           Tail_ESS
sigma_Intercept                2885
alpha_Intercept                3568
Intercept                      2387
patient1                       2257
condN                          1890
armND                          1852
patient1:condN                 3400
patient1:armND                 2923
condN:armND                    3083
patient1:condN:armND           3163
sigma_patient1                 3217
sigma_condN                    3166
sigma_armND                    2540
sigma_patient1:condN           2930
sigma_patient1:armND           2991
sigma_condN:armND              2775
sigma_patient1:condN:armND     2800
alpha_patient1                 2958
alpha_condN                    3268
alpha_armND                    3073
alpha_patient1:condN           2985
alpha_patient1:armND           3109
alpha_condN:armND              2960
alpha_patient1:condN:armND     2845

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(m11_transFlat, ndraws = 100)

Fig. 31. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "mu"),
                      plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "sigma"),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "alpha"),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transFlat,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 32. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.


2.2.2 Conclusions

wrap_plots(plot(conditional_effects(m11_transFlat,
                                    dpar = "mu",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]] + 
             scale_color_colorblind() +
             theme(legend.position = "none"),
           plot(conditional_effects(m11_transFlat,
                                    dpar = "sigma",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]] +
             scale_color_colorblind()  + 
             theme(legend.position = "none"), 
           plot(conditional_effects(m11_transFlat,
                                    dpar = "alpha",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]]) +
  scale_color_colorblind() + 
  plot_annotation(tag_levels = "A")

Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the full model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).


In transFlat there is a small but most likely positive difference between controls and patients, but all the other main effects are not clearly different from zero (Fig. 30A).

Variation shows a curiously extreme estimated value for healthy subjects playing controlled strokes on the non-dominant side (Fig. 30B) which seems to be mostly driven by a relatively large estimated value for the non-dominant side (sigma_armND in Tab. 14).

Fig. 30C shows that the left tail of transFlat’s distribution gets heavier in the controlled condition, suggesting a tendency in players to produce a larger amount of low-valued transFlat values when performing unfamiliar strokes.

Since there is no clear difference in the estimated location parameter \(\mu\) for transFlat values between conditions, this descriptor is probably not a suitable cue to distinguish between conditions in a listening test.

2.3 Univariate Models – transCrest

2.3.1 The Full Model: Modeling the Interactions of Group, Instruction, and Side on Shape and Spread

(m11_transCrest_form <-bf(transCrest ~ 0 + Intercept + patient * cond * arm + 
                            (1 | subj) + 
                            (1 | cond) +
                            (1 | arm),
                          sigma ~ patient * cond * arm,
                          alpha ~ patient * cond * arm))
transCrest ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
sigma ~ patient * cond * arm
alpha ~ patient * cond * arm
(m11_transCrest_prior <- c(set_prior("normal(2.5, 10)",    class = "b", coef = "Intercept"),
                           set_prior("normal(0, 10)",  class = "b"),
                           set_prior("normal(0, 5)", class = "b", coef = "armND"),
                           set_prior("normal(0, 5)", class = "b", coef = "condN"),
                           set_prior("normal(0, 0.01)",      class = "sd"),
                           set_prior("normal(0, 5)",      class = "Intercept", dpar = "sigma"),
                           set_prior("normal(0, 5)",      class = "b",         dpar = "sigma"),
                           set_prior("normal(0, 2)",      class = "Intercept", dpar = "alpha"),
                           set_prior("normal(0, 2)",      class = "b",         dpar = "alpha")
))
           prior     class      coef group resp  dpar nlpar   lb   ub source
 normal(2.5, 10)         b Intercept                        <NA> <NA>   user
   normal(0, 10)         b                                  <NA> <NA>   user
    normal(0, 5)         b     armND                        <NA> <NA>   user
    normal(0, 5)         b     condN                        <NA> <NA>   user
 normal(0, 0.01)        sd                                  <NA> <NA>   user
    normal(0, 5) Intercept                      sigma       <NA> <NA>   user
    normal(0, 5)         b                      sigma       <NA> <NA>   user
    normal(0, 2) Intercept                      alpha       <NA> <NA>   user
    normal(0, 2)         b                      alpha       <NA> <NA>   user
if (MODEL) {
  m11_transCrest <- brm(m11_transCrest_form,
                        prior = m11_transCrest_prior,
                        family = skew_normal(),
                        init = "0",
                        control = list(max_treedepth = 15),
                        data = drum_beats)
  m11_transCrest <- add_criterion(m11_transCrest, 
                                  "loo",
                                  reloo = TRUE)
  save(m11_transCrest, 
       file = "m11_transCrest.rda")
} else {
  load("m11_transCrest.rda")
}

Tab. 15. Model summary.

(m11_transCrest_summary <- summary(m11_transCrest, 
                                   priors = TRUE))
 Family: skew_normal 
  Links: mu = identity; sigma = log; alpha = identity 
Formula: transCrest ~ 0 + Intercept + patient * cond * arm + (1 | subj) + (1 | cond) + (1 | arm) 
         sigma ~ patient * cond * arm
         alpha ~ patient * cond * arm
   Data: drum_beats (Number of observations: 1102) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Priors:
b ~ normal(0, 10)
b_armND ~ normal(0, 5)
b_condN ~ normal(0, 5)
b_Intercept ~ normal(2.5, 10)
b_alpha ~ normal(0, 2)
b_sigma ~ normal(0, 5)
Intercept_alpha ~ normal(0, 2)
Intercept_sigma ~ normal(0, 5)
<lower=0> sd ~ normal(0, 0.01)

Multilevel Hyperparameters:
~arm (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.01      0.01     0.00     0.02 1.00     3013     1546

~cond (Number of levels: 2) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.01      0.01     0.00     0.02 1.00     3522     1983

~subj (Number of levels: 8) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.09      0.01     0.08     0.10 1.00     4434     2804

Regression Coefficients:
                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sigma_Intercept               -0.96      0.06    -1.08    -0.83 1.00     2792
alpha_Intercept                5.57      0.91     3.94     7.53 1.00     3944
Intercept                      3.03      0.06     2.92     3.15 1.00     2042
patient1                      -0.51      0.08    -0.66    -0.36 1.00     1965
condN                         -0.15      0.04    -0.23    -0.08 1.00     2481
armND                         -0.18      0.05    -0.27    -0.09 1.00     2511
patient1:condN                 0.08      0.05    -0.01     0.16 1.00     2428
patient1:armND                 0.30      0.06     0.18     0.41 1.00     2313
condN:armND                    0.10      0.06    -0.02     0.22 1.00     2326
patient1:condN:armND          -0.18      0.07    -0.33    -0.04 1.00     2269
sigma_patient1                -0.43      0.10    -0.62    -0.24 1.00     2580
sigma_condN                   -0.63      0.09    -0.80    -0.45 1.00     2853
sigma_armND                    0.05      0.09    -0.12     0.22 1.00     2428
sigma_patient1:condN           0.49      0.14     0.21     0.75 1.00     2631
sigma_patient1:armND           0.42      0.14     0.15     0.69 1.00     2213
sigma_condN:armND              0.63      0.12     0.39     0.87 1.00     2378
sigma_patient1:condN:armND    -0.91      0.18    -1.26    -0.54 1.00     2342
alpha_patient1                -0.06      1.17    -2.32     2.24 1.00     3764
alpha_condN                   -3.19      0.98    -5.15    -1.29 1.00     3858
alpha_armND                    1.04      1.26    -1.43     3.53 1.00     3746
alpha_patient1:condN           0.54      1.19    -1.87     2.88 1.00     3727
alpha_patient1:armND           3.86      1.51     0.98     6.80 1.00     4501
alpha_condN:armND              0.22      1.34    -2.37     2.92 1.00     3836
alpha_patient1:condN:armND     1.51      1.71    -1.75     4.86 1.00     5090
                           Tail_ESS
sigma_Intercept                3223
alpha_Intercept                3185
Intercept                      2562
patient1                       2268
condN                          2898
armND                          2964
patient1:condN                 2827
patient1:armND                 2950
condN:armND                    2956
patient1:condN:armND           2671
sigma_patient1                 3033
sigma_condN                    3190
sigma_armND                    2700
sigma_patient1:condN           2689
sigma_patient1:armND           2574
sigma_condN:armND              3004
sigma_patient1:condN:armND     2649
alpha_patient1                 2774
alpha_condN                    3036
alpha_armND                    2685
alpha_patient1:condN           2809
alpha_patient1:armND           2981
alpha_condN:armND              2799
alpha_patient1:condN:armND     3189

Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pp_check(m11_transCrest, ndraws = 100)

Fig. 34. Posterior predictive check. The thick blue line shows the distribution of the empirical data. The thin blue lines are one-hundred realizations of data generated from parameters estimated by the model.


beautify_my_plot(plot(conditional_effects(m11_transCrest,
                    dpar = "mu"),
                    plot = FALSE))

conditions <- make_conditions(drum_beats, 
                              vars = "arm")
beautify_my_plot(plot(conditional_effects(m11_transCrest,
                                          dpar = "mu",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transCrest,
                    dpar = "sigma"),
                    plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transCrest,
                                          dpar = "sigma",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transCrest,
                                          dpar = "alpha"),
                      plot = FALSE))

beautify_my_plot(plot(conditional_effects(m11_transCrest,
                                          dpar = "alpha",
                                          effects = "patient:cond",
                                          conditions = conditions),
                      plot = FALSE))

Fig. 35. Conditional plots. Had to take out all the fancy figure polishing, subplot lettering etc. Subplots are supposed to be numbered from A to … starting in the top left corner, proceeding row-wise.


2.3.2 Conclusions

wrap_plots(plot(conditional_effects(m11_transCrest,
                                    dpar = "mu",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]] + 
             scale_color_colorblind() +
             theme(legend.position = "none"),
           plot(conditional_effects(m11_transCrest,
                                    dpar = "sigma",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]] +
             scale_color_colorblind()  + 
             theme(legend.position = "none"), 
           plot(conditional_effects(m11_transCrest,
                                    dpar = "alpha",
                                    effects = "patient:cond",
                                    conditions = conditions),
                plot = FALSE)[[1]]) +
  scale_color_colorblind() + 
  plot_annotation(tag_levels = "A")

Fig. 30. Three-way interaction plots of the three modeled parameters \(\mu\), \(\sigma\), and \(\alpha\) as estimated by the best-fitting model. The x axis shows group membership (0 = healthy subjects, 1 = patients). The instruction conditions are color-coded (black = controlled strokes, brown = normal strokes), and the facets show the side (1 = dominant, 2 = non-dominant).


The location parameter \(\mu\) of the transCrest data has been estimated to be clearly different between groups, conditions, and sides such that there is even a three-way interaction between the main effects (Fig. 30A): patients have low \(\mu\) on the dominant side and higher \(\mu\) on the non-dominant side, while the opposite holds for healthy controls.

The estimation of spread displays some weird patterns, with \(\sigam\) generally being larger in the in the controlled instruction condition, but with one strangely small estimate in patients on the dominant side. The opposite holds for the normal instruction condition, with a saliently large value in healthy subjects when playing normal strokes with the non-dominant arm.

The distributio of the data has a heavy tail on the right side seems to follow a pattern of “normal” < “controlled”, and “dominant” < “non-dominat”. Only for patients on the non-dominant side skewness seems to be about the same, irrespective of instruction.

With the three-way interaction in the location parameter this descriptor seems a promising candidate for a listening test because if we aim for a single descriptor to be sufficient for listeners to reliably distinguish between all three factors, only a descriptor satisfying this interaction can be used to delineate the different strokes. But since the original aim–as I understand it–is to primarily distinguish between controlled and normal strokes, every descriptor with a clear difference between these two factor levels would be a good candidate from a statistical point of view. In the three descriptors investigated here, transSC most clearly distinguishes (again: from a stats point of view) between normal and controlled instruction conditions.

3 System Setup

Data wrangling and analyses were carried out with the statistical package R [R version 4.4.2 (2024-10-31); R Core Team (2020)]. Bayesian modeling was done with the package brms (Bürkner 2017, 2018) which uses the probabilistic language Stan as back end (Carpenter et al. 2017). Plots were done with the packages bayesplot (Gabry and Mahr 2020) and ggplot2 (Wickham 2016).

References

Bates, D. M. 2010. Lme4: Mixed-Effects Modeling with r. Springer. http://lme4.r-forge.r-project.org/lMMwR/lrgprt.pdf.
Bürkner, Paul-Christian. 2017. “Brms: An r Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1). https://doi.org/10.18637/jss.v080.i01.
———. 2018. “Advanced Bayesian Multilevel Modeling with the R Package Brms.” The R Journal 10 (1).
Câmara, Guilherme Schmidt, Kristian Nymoen, Olivier Lartillot, and Anne Danielsen. 2020. “Effects of Instructed Timing on Electric Guitar and Bass Sound in Groove Performance.” The Journal of the Acoustical Society of America 147 (2): 1028–41. https://doi.org/10.1121/10.0000724.
Carpenter, Bob, Andrew Gelman, Matthew Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “Stan: A Probabilistic Programming Language.” Journal of Statistical Software, Articles 76 (1): 1–32. https://doi.org/10.18637/jss.v076.i01.
Danielsen, Anne, Carl Haakon Waadeland, Henrik G. Sundt, and Maria A. G. Witek. 2015. “Effects of Instructed Timing and Tempo on Snare Drum Sound in Drum Kit Performance.” The Journal of the Acoustical Society of America 138 (4): 2301–16. https://doi.org/10.1121/1.4930950.
Gabry, Jonah, and Tristan Mahr. 2020. “Bayesplot: Plotting for Bayesian Models.” https://mc-stan.org/bayesplot.
R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Tibshirani, Robert. 1996. “Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society: Series B (Methodological) 58 (1): 267–88. https://doi.org/10.1111/j.2517-6161.1996.tb02080.x.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.
Wilkinson, G. N., and C. E. Rogers. 1973. “Symbolic Description of Factorial Models for Analysis of Variance.” Applied Statistics 22 (3): 392–99. https://doi.org/10.2307/2346786.